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We develop a theoretical description of intravalley scattering of quasiparticles in graphene from multiple 
short-range scatterers of size much greater than the carbon-carbon bond length. Our theory provides a method to 
rapidly calculate the Green's function in graphene for arbitrary configurations of scatterers. We demonstrate that 
non-collinear multiple scattering trajectories generate pseudospin rotations that alter quasiparticle interference, 
resulting in significant modifications to the shape, intensity, and pattern of the interference fringes in the local 
density of states (LDOS). We illustrate these effects via theoretical calculations of the LDOS for a variety of 
scattering configurations in single layer graphene. A clear understanding of impurity scattering in graphene is a 
step towards exploiting graphene's unique properties to build future devices. 

PACS numbers: 



I. INTRODUCTION 

Due to the fact that the graphene unit cell in single layer graphene consists of two inequivalent carbon atoms [the A and B car- 
bon atoms shown in Fig. [TjA)], the quasiparticles in graphene can be described as a pseudospin. This imbues the quasiparticles 
in graphene with a sense of chirality that depends upon the pseudospin direction. Compared to conventional two-dimensional 
electron gases (2DEGs), graphene's chiral nature provides it with tremendous potential for use in ultrafast electronic and pseu- 
dospintronic devices^^ An important question that has been extensively studied is the effect of point defects or point scatterers 
upon the local density of states (LDOS) in graphene, where the size of the "point" scatterer is on the order of or smaller than 
the C — C bond, b = 0.142 nm^ ^. Scattering from defects of this size gives rise to both intravalley and intervalley scattering, 
the latter of which significantly affects the transport properties in graphene. Most studies have focused on calculating the effects 
of scattering from one or two point defects (a quantum corral of point defects in graphene has also been studied^) on the LDOS 
and have used the scattering amplitudes from single point defects as input in modeling the mobilities in graphene^^. In these 
models, the calculated LDOS is sensitive to the particular placement of the scatterer(s) on the graphene lattice. However, as 
the size or scattering length of the scatterers becomes much larger than the C — C bond length, the LDOS should become less 
sensitive to the actual placement of the scatterers on the graphene lattice and should depend only on the relative configuration of 
scatterers especially at low wavelengths ^ <C ^ . While controlling the distribution of point defects in a device is challenging, 
the controlled introduction and arrangement of large scatterers in graphene shoul d be e asier to accomplish experimentally by, for 
instance, placing small metallic islands placed atop or below a graphene surfacd^^^. In this case, such large scatterers should 
predominately cause intravalley scattering. 

In this work, we present a theory for intravalley scattering of quasiparticles from multiple scatterers when the wavelength of 
the quasiparticles is comparable to the size of the scatterers. The ^— matrix operators for all partial waves for a single scatterer are 
derived and used in developing a theory of multiple scattering in graphene. Our calculations suggest many partial waves must be 
included in order to accurately calculate the LDOS even when the quasiparticle 's wavelength exceeds the scattering length. That 
is, the usual ^-wave approximation, which is ubiquitous in low-energy scattering of free particles, breaks down in graphene. We 
calculate both the Green's function and the LDOS in graphene in the presence of multiple scatterers. From our calculations, the 
chiral or pseudospin nature of graphene results in changes to the intensity, shape, and pattern of quasiparticle interference in the 
LDOS compared to the LDOS found in a 2DEG. Our calculations should be experimentally verifiable using scanning tunneling 
microscopy (STM) or in "simulations" of graphene systems, such as in the scattering of an atom from multiple atoms confined 
in a hexagonal optical latticJ^. 

II. THEORY 
A. Graphene preliminaries 

In this section, we briefly review the physics of graphene^ that will be relevant to the scattering theory subsequently developed 
in this paper. A single layer of graphene consists of two displaced triangular lattices of carbon atoms, A and B, that generate a 
hexagonal or honeycomb structure as shown in Fig.[TjA). Each unit cell of graphene consists of an A and B lattice site where 
the unit cell is defined by two integers, j = [m^n], and the positions of the A and B lattice sites in the unit cell are 

fj = ma^ + na- and = rj — by respectively, where a± = ± ^^x-\- ^y are the lattice vectors for the honeycomb lattice and 
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Z? = 1.42 A is the carbon-carbon bond length. The corresponding reciprocal lattice vectors are given by = || (±>/3£+j), 
with a'^-a±= 271 and a'^-a^ = 0. 

Defining the following lattice wave functions over Niat unit celliEl 



(1) 



where <^>a(b);(|'^ — '^^^^D denotes an orbital centered on the A(B) lattice site in the unit cell. Tight-binding calculation^ 
taking into account nearest neighbor coupling have previously demonstrated that there exist four distinct, zero-energy states in 



graphene that are given by 



^ ((^f{r,K)±^l\r,K)^ and ^ (^^^^\r,-K)±^f{j,-K)^, where ±K = i?^^ = ±=^£ 

The dispersion relation for graphene is linear when expanded about ±K for small k [\k\ ^ |^|]. In this case, the lattice wave 
functions aik±K can be written in terms of the lattice wave functions at <l>^^[^) (r, ±^ + ^) : 



^fir, K) ± ^fir, K) ] and ^ {^f{r, -K) ± ^f{r, -K)^ , where ±K = ±- 



e"'-'^f^^^{f,±K). The ener- 
gies and eigenstates of the tight-binding Hamiltonian about ±K can be written as E^^ 

c^{f,k)<$>B{^^ ^K) respectively, where the coefficients or "envelope" functions, c^^^^(r,X), and energies are determined by solv- 
ing the following equations: 



d_ 
dr 



dr 



i_d_ 

~rde 



= ±hVFkL. 



(2) 



where hvp = 1.0558 x 10~^^J-m, and L± = jj^e^'^ (f ^ r ^) = ^ ^^|;)' solutions to Eq. (2) are parameterized by 
the wave vector k = {kcos O^^ksinG^) with |^| <C |^| and are given by 



1 
1 

7=2 



Jk-r 



= e 



jk-r I 



Jk-r 



^(<DA(r,^)±e'^*<I>B(r,^) 



Jk-r 



±)_^ = ^ (<Da(?, -K) ± e'^*<I>B(r, -K) 



(3) 



where the energies are Ej^ = ±hvFk for ^^^^^(r,^,±) and E_^ = ^fivph for ^^i|^(r,^,±) [the superscript "clean" refers to 

graphene in the absence of scatterers] . The linear dispersion relation about gives the same results as the exact tight-binding 
calculations for \k\b < 0.2. Such a linear dispersion relation is analogous to the dispersion relation found for massless Dirac 
fermions^^ and therefore the wave vectors, are often referred to as Dirac points. 

From Eq. (pj), the envelope functions for ^^^^ are plane waves. Since we are developing a theory for scattering from localized 
scatterers, we can express the eigenstates for the graphene Hamiltonian in cylindrical coordinates, which are given by: 



MO 



Hl''^\kr)e"'^ 
±iHl\:f\kr)e'('+'^<^ 



(4) 



J K 



about the K Dirac point [the cylindrical eigenstates about the —K Dirac point are found by simply replacing K by —K in Eq. (Wl], 

mid Hp^\kr) are Hankel functic 
waves about r = 0, respectively. 



and Hj^'^\kr) are Hankel functions of order /. The states J^^^^-^e^^^ and ^j^^\e^^^ represent outgoing and incoming cylindrical 



B. The matrix for intravalley scattering 



In the following, a theory for intravalley scattering 



k±K^k'±K 



from cylindrically symmetric scatterers placed atop a 



graphene sheet will be developed. First, we consider the case of a single, cylindrically symmetrical scatterer of radius a centered 
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FIG. 1: [Color Online] (A) Honeycomb lattice structure of monolayer graphene consisting of two triangular lattices [denoted by A (red) and 
B (blue) carbon atoms] that are shifted relative to each other by —by where b = \ A2 A is the C — C bond length. Each unit cell consists 
of an A and B lattice site [a typical unit cell is denoted by the dotted ellipse in Fig. [iJa)], with the position of the unit cell given by 

fj = md^ -\-nd- where j = [m,n] and d± = ±^^£+ are the lattice vectors. (B) Scattering configuration in graphene for an incoming 
"wave", <l>^^(r,X, ±) = ^^^'^\^)±k incident upon a scatterer located at r„. 



at rn and represented by a potential acting equally on both the A and B lattice site^, V{r) = Vq for |r — r„| <a and V{f) = 
for \f — fn\ > a; that is, the scatterer behaves as a uniform potential of radius |r — r„| < a. Scattering from such potentials in 



graphene has been previously studiecP^^ISl Pqj. intervalley scattering processes 



k±K 



<C Tla^, which implies that the radius of the scatterer should satisfy 



2n 



to be neglected, we require 
3.7 A. 



We will first consider intravalley scattering of a free particle wave with E >Q about each Dirac point, <J>^^(r,^, ±) = 

^^^^(r,^, ±), incident upon a scatterer located at r„ [Fig. I|b)]. In this case, the full scattered wavefunction can be written as 
follows: 



where = |r — r„|, e^^^"" = 
scattered wave function is given by: 



{r-rn)-x-^i{r-rn)-y 



Pn 



, and 5/ are phase shifts of the outgoing cylindrical partial waves, ^^^L The 



(5) 



where si = ^ ' ^ ^ the scattering amplitude of the /^^ partial wave, which is determined by continuity of the wavefunction at 
pn = ci and is given by 



Ji{k' a)JiJ^i{ka) - Ji(ka)JiJ^i(k' a) 



(k'a)H\^^ {ka) - Ji (k'a)H\'^^ (ka) 



(1) 



(6) 



where k' = k- 



^ and // is a bessel function of order /. Note that si is the same for scattering about K and —K and satisfies the 



unitarity condition — Re[^/] = |^/| for all /. 

Using the above results, we now derive the ^— matrix for scattering of arbitrary incident waves of energy E from a scat- 
terer at Vyi. The derivation follows that used for deriving the ^— matrix in two-dimensional electron gases with Rashba 
spin— orbit coupling^^. Writing the incident wave at the scatterer as <I>^^(r„,^, ±) = ^^^'^^ and using the fact that 



±^(%'«l%'i5)±^ = ^ap ^^^^sgniif ^ '*^here sgn(/) gives the sig: 



can be rewritten as: 



e'''"J^^%{Pn:k,±hi{9^,,±\\\ej^,±)^g 



= I 



;l c.,JW„ 



±iH\'_^^{kpn)e'^'' iHll\{kpn)e' 





'i'e-'V-^"' 


' ±K 





I 



Hj'\kpn) ±Hf\kp„)L. 



2 \^ ±iHl^,{kpn)e^''' iH\%{kpn)e'^''L. ^ 



oo ./ 

E - 



(1), 



2 \ ±/H« (;:p„)e'^'' H/''(^:p„) 



±K 



y^^"^^gn(-') 



,^iL_n*L"^(r„,l±) (7) 



Eq. Q can be regrouped into terms with the same scattering ampHtude, si [note from Eq. ([6]), si = Using the fact that 

any wavefunction, that is a solution to Eq. ([2]) satisfies the following relations: 



1 tl. 



^ - = 



Eq. (|7]) can be rewritten as: 



/=0 



where 



/ Hj^\kpn)e'^^'- ±iHll\{kpn)e-'^'+^^^- 
'4hVF \±iHjl\{kp„)e'^'+'^^^- Hj^\kp„)e-''^" 

_ik_ ( V^[Hl^\kpn)] TiLL[H^ll(kpn)e-'^''] ' 
1^ ±/L'+[//r^(^P„)e''"] L'_[//^''(^p„)] , 



±K 



±K 



and 7] is the /-partial wave f— matrix operator given by: 



(8) 



(9) 



(10) 



T - — 



IJ 

lL 



±K 



(11) 



Note that Go^±K{r^fn^E) is simply the Green's function for single layer graphene about the ±K Dirac poinP^. 

In a similar manner, the total Green's function, G^^{f^f' ,E), can be calculated and is given as [using the above notation]: 



\G2i{r,r',E) G22{r,r' ,E) 
Aihvpsi 



±K 



= ^o,±l(^'^''^) + L — I — ^/,±^(^'^"'^)^/,±l G^^(rn,r,E) 



1=0 



1=0 



For a single scatterer at r„, 

LL[Gn{rr,,r,E)] LL[Gn{rn.r' ,E)] 

rl 



i-iy+^k ( Hl^\kp'„)e-'"^!' TiHll\{kp'„)e-'^'+^^'^!' 



4hVF \^iHY^,{kp',)e 



H\^\kp'„)e'^^'-' 



(12) 



(13) 



±K 
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where = |r„ — r'| and e^^^^^ = ^ 



{r'-rny{x±iy) 



P'n 



. In Eq. (13 



we used the relation that L 



Inserting Eq. ([13]) into Eq. ([12]) completely determines the Green's function in the presence of a single scatterer. 

The LDOS, p^^{f^E), which is an important quantity that has been previously measured in STM experiments on 

grapheneSEU^ calculated from the Green's function using the relation^^ p^^{f^E) = — ^Im 



pj^^ir^E) for single layer graphene in the presence of a scatterer at r„ is given by: 
k 



G^iir,r,E) 



. With Eq. (12), 




{H\'\kpn)y 

{H\'\kPn)y 



{Hil\{kPn)) 
{Hil\{kPn)) 




±K 



pff{r,E) + 5p^^{r,E) 



■J clean/' 
±K 



p^_^Y\r,E) 



(14) 



where ^Jf^(r, ±:K) and ^^^(r, ^K) are the lattice wave functions [Eq. ([l|)] evaluated at r, 



pff "(r,£) - 



is the LDOS for single layer graphene about the ±:K Dirac points in the 



absence of any scatterers, 5p^^(r,£) is the change in the LDOS due to the scatterer, and 



5p^^{r,E) 
ptT^r.E) 



= £lm L {^{H['\kpn))\{H['^,{kpn^^ 



1=0 



(15) 



is the "envelope" function of the Friedel oscillations in the LDOS. 



For/:p,>l, (Hl^\kpn)) - 
5p^^{r,E) 



2 ^2ikpn —i^—iliz _j_ Q 



and so in this limit, Eq. ( 15 ) becomes: 



Elm 



/ / ilkpn _ J2kpn 



o 



1 

kpn 



= 




(16) 



lations in 2D electron gas (2DEG) or achiral systemPSEll decay as ^Lo^m ^isi (^H\^\kpn)^ ^ l4=o^^ [si{-\)^ e^'^P^] . 

[ Fig.j2|A). From Eq. ( 16), the {kpn)~^ dependence 



Thus the Friedel oscillations in Sp^^ decay as (kpn) ^ in graphene^, whereas for kpn 1, the envelope of the Friedel oscil- 



This feature is illustrated in the numerical simulation of h£1?5L shown in 

in the Friedel oscillations is due to the destructive interference between H^^\kpn) and H^^^{kpn) along a collinear scattering 
trajectory from the position r, to the scatterer and then back [r ^ r„ ^ r]. While it is tempting to attribute the {kpn)~^ depen- 
dence of the Friedel oscillations in graphene to the absenc^of intravalley backscattering (k —k), this can not be the complete 
explanation behind the decay found in Eq. ([16]). As a case in point, the calculated LDOS [using the multiple scattering theory 
developed in the next section] outside an elliptical array of = 30 scatterers [semimajor and semiminor axis of 20 nm and 10 
nm respectively] and along the direction of the semimajor axis is shown in Fig.[2[B). Although backscattering (k —k] from 
an elliptical array of scatterers is still prohibited in graphene by time-reversal symmetr}^, the Friedel oscillations for graphene 
extend out a hundred nanometers from the scatterer unlike that found for a single scatterer [Fig. ^A)] where the Friedel os- 
cillations would have died out within tens of nanometers from the scatterer array. For the elliptical array of scatterers, there 

exist noncollinear multiple scattering trajectories that prevent the complete destructive interference between the H^^\kpn) and 
(kpn) components of the /^^ -partial wave, thereby leading to a slower decay of the Friedel oscillations. 



III. MULTIPLE SCATTERING IN GRAPHENE 



The above theory can be extended to investigate multiple scattering from N scatterers in graphene. However, for multiple 
scatterers, it is computationally unfeasible to consider all partial waves. Fortunately for most physical potentials, the scattering 
amplitude, su is nonneglible for only a small subset of /, thereby justifying the use of only a few partial waves in the scattering 
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|r-rn|(nm) 



FIG. 2: [Color Online] The calculated relative change in the LDOS, dean A L in the presence of (A) a single scatterer at r„ and (B) an elliptical 

P±k ^^'^^ 

arrangement ofN = 30 identical scatterers [major and minor axes of 20 nm and 10 nm respectively] in both graphene [solid, blue curve] and 
in a 2D electron gas (2DEG) or achiral system [dotted, red curve] as a function of (A) the distance from the single scatter or as a function of 
(B) the distance from the outermost scatterer along the semimajor axis of the ellipse. In both calculations, a = 10 A, Vq = 4 eV, the scattering 
amplitudes were given by si in Eq. ([gJ, and k = 0.485 nm~^ [corresponding to E = 0.32 eV for graphene]. In both (A) and (B), the first five 
partial waves [/ = to / = 4] were included in the calculations, which provides an accuracy of the numerical results to within 1 percent. (A) 
For the single scatterer case, the Friedel oscillations in the LDOS decrease as (kpn)~^ for the graphene/chiral case [blue, solid curve], whereas 
they decrease as {kpn)~^ in a 2DEG/achiral case [red, dotted curve]. (B) Due to multiple scattering within the elliptical array of scatterers, the 
Friedel oscillations in the LDOS decrease at roughly the same rate as those found in a 2DEG and can be observed at distances up to 100 nm 
from the scatterers. 



calculations. In the following theory, only the first lynax + 1 partial waves [/ = to / = lmax\ will be considered [a discussion of 
the proper choice of Imax will be provided later in the paper] . Denoting the position of the scatterer by fj and the partial 

wave scattering amplitude from scatterer j as s\^\ the total Green's function at energy E = hVpk > 0, can be written as [r, 7^ rj 
for 7 = 1 to j = N]\ 



U) 



G^^{ry,E) = Go ,^(r,r',^) + ^^ -^Gj^^{r,rj,E)Tj^^ G^^{rj,r' ,E) 



j=u=o 



(17) 



In Eq. ( 17 ), knowledge of 7] Gj^^ijj^r'^E) for each partial wave / = to / = Imax and for each scatterer 7 = 1 to j = N 
completely determines the total Green's function, G^^{f^r',E). In this case, there are 4N{l}nax-^ 1) unknowns that can be 



determined self-consistently from the corresponding Foldy— La^J^SET equations derived from Eq. (17 ): 



Imax /^ifi^ps 



U) 



j^nl=0 



fovn = I to n = N and /' = to /' = Lax- where 



(18) 



ie„j{i+i'+i) 



±K 



Defining the following 2{l,nax + 1) x 2 matrix for each scatterer j: 



TG (o-,/)^^ = 



/ 



V \a.,±K[Go^±i^7j,r' ,E)] j 



(19) 



(20) 



and the following 2N{ljnax + 1) x 2 matrices: 



/ rG±^(ri,r') \ 



/ rGo,±^(ri,r') \ 



Eqs. ( 18 ( can be written compactly for n = 1 to n = and for /' = to /' = Imax as: 



TG^^(/) = (1-TT) TGo^±^(r') 



where 1 is the IN [Imax + 1) x IN [Imax + 1) identity matrix and f T is a 2N{lmax + 1) x IN {Imax + 1) matrix given by: 



TT 



\TT±k(^nA) TT^^ijN.ri) TT^^ijN.h) 



where 



(21) 



(22) 



(23) 









... ^ 












7 - 




7 - 












7 - 


^2,±^(^^' 0') 


7 - 
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7 


^2,±^(^"'0") 


... 7 

Imaxi^K 


_^lmax.±Ki^'''^j\ 



and 



( s\j^ 

^[-^"^ 



0^: 



\ 






ii) 



V 

Finally, Gj^^ir^r' ^E) in Eq. |l7| ) can be written compactly as: 

G^^ir^y^E) = Go_^^(r,r',£)+GG(r)fG(r') 



where GG(r) is a 2 x 2N{l,„ax + 1) matrix, GG(r 



1 
1 



GG{r,rj) 



Go^±^(r,r',£) + GG(r)(^l-TTj TGo(/) 
GG(r, ri ) GG(r, r2) ... GG(r, Fat) where 



(25) 



(26) 



(27) 



Using G_|_^ (?,?',£") in Eq. (26), the LDOS, p_^_^{r,E), and all transport properties can be calculated in graphene in the presence 
of N scatterers. 
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IV. RESULTS AND DISCUSSION 



In this section, we use Eq. (26) to calculate the effects of quasiparticle interference on the LDOS in single layer graphene 
in the presence of multiple scatterers. As the seminal work of the Eigler group^^ demonstrated, quantum corrals comprised of 
adatoms placed atop a metal surface lead to quasiparticle interference, where changes in the LDOS mimic the behavior of an 
electron confined in a billiard with leaky walls^^"^. While the quasiparticles in graphene cannot be confined by such corrals 
due to their Dirac-like nature (this phenomenon is a consequence of the Klein paradox^^), non-collinear multiple scattering 
trajectories in graphene can result in dramatic changes in the LDOS compared with single scattering trajectories, which was 
previouslyillustrated in Fig.|2] 

In Fig. 3 the relative change in the LDOS, fifai^^L , is shown for a rectangular array [10 nm by 20 nm] of A/^ = 36 identical 



scatterersla = 1 nm and V = 4 eV] in both graphene [left, calculated using Eq. ([26])] and in a 2DEG [right, calculated from 
Eq. ( |A8| )]. For better comparison between the 2DEG and graphene cases, the same magnitude of wave vector, k, and scattering 
amplitudes, si [Eq. ([6])] were used in evaluating the LDOS for both the graphene [Eq. (26)] and the 2DEG [Eq. (AS)] cases. 
Therefore, the differences between calculated LDOS for the graphene and 2DEG cases arise solely from differences in the 
"pure" Green's function in the absence of scatterers, the chiral Gq ^^(r, r',£) ^ G^^ ^{f' ^r^E) in Eq. ( |iq| ) for graphene vs. the 

achiral Go{f^f'^E) = Go{f',f^E) = ^^^^|r — r'|^ for the 2DEG. For accurate calculations of Gir^f^E) and hence 

the LDOS^ r, all scattering trajectories beginning and ending at r must be taken into account when closely spaced scatterers 
are present^^, as is done in Eq. ( |26| ) and Eq. ( A8 ). 

From Fig. [3] a few general observations can be made. First, the Friedel oscillations outside of the scatterer array are in general 
weaker in graphene than those found in the 2DEG, which is due to the fact that collinear backscattering trajectories [Fig.|2jA)] 
contribute little to the LDOS away from the scatterers. Second, the shape of the Friedel oscillations outside the rectangular array 
of scatterers differs significantly between graphene and the 2DEG. For the 2DEG, the Friedel oscillations are rectangular in shape 
away from the array of scatterers. In graphene, however, the shape of the Friedel oscillations is not rectangular and is highly 
anisotropic. Treating <I>^^(r, and <I>^^(r, as a pseudospin pair, multiple scattering trajectories result in pseudospin 
rotations that alter the quasiparticle intereference and thereby modulate the LDOS since the interference between two different 
trajectories is maximal only when the pseudospins from each trajectory are pointed along the same direction. A similar result was 
predicted for the LDOS of quantum corrals with Rashba spin— orbit coupling^^, where multiple scattering trajectories generate 
noncommuting rotations that result in extra modulations in the LDOS compared to the LDOS found in 2DEGs. Lastly, the LDOS 
inside the rectangular array is also quite different between graphene and the 2DEG case. While previous work on a 2DEG has 
demonstrated that the change in the LDOS can be quite large inside a quantum corral due to quasiparticle interference^^, these 
confinement resonances are, in general, not observed for scattering in graphene corrals, which is a consequence of the inability of 
the scattering potential to confine the quasiparticles in graphene [from calculations, resonances were not observed in a circular 
array of ^-wave scatterers, but some resonances were observed as more partial waves were included in the calculation^^]. It 
should be noted that since the scatterers used in the above theory affect both the A and B lattice sites identically, the LDOS 
calculated from Eq. (26) and shown in Fig. |3] would be valid description of low-resolution/coarse-grained STM images with 
spatial resolutions larger than than the C — C bond length, b= 1.42 A. 

In the calculations shown in Fig. 3 l^ax = 4 was chosen in order to achieve a relative accuracy in the fifai?5L profiles 

U _ P±K ^"'^^ 

of one percent, even though, on average, |^o| ^ \^i^o\ [in Fig.pl |^o| ^ 0.94, |^i| ^ L5 x 10 ^I^s-qI, |^2| ^ 6.5 x 10 ^I^-qI, 
\s3\^3 X 10~^, and |^4| ^ 6 x 10~^|^o|]- However, what matters in the calculation of G^^{f^f^ ^E) in Eq. (26) is the inverse of 

XT, which contains off-diagonal terms that are proportional to H^^\krnm)si' that can be quite large when krnm ^ 1- Therefore, 
as the distance between scatterers decreases, more partial waves must be included in the calculation. This is due to the fact that 
when the distance between scatterers is comparable to their scattering length/size, more partial waves are needed in order to 
insure that the total wave function satisfies the boundary conditions at the edge of each scatterer [p„ = a]. In Fig. [3] the nearest 
neighbor distance between scatterers was 33.33 A^1.6x2a, which required Imax = 4 to achieve a relative accuracy of one 
percent. This is illlustrated in Fig. |4jC), where the calculations of the LDOS for an elliptical array [semimajor and semiminor 
axis of 20 nm and 10 nm respectively] of A/^ = 30 scatterers [a= I nm and V = 4 eV] are shown for Imax = 4 [left] and Imax = 
[right]. In this calculation, k = 0.424nm~^ and |^o| = 0.9319 |^i | = 0.0147; however, the inclusion of the first l^ax + 1=5 
partial waves results in small changes to the LDOS that can be seen around the elliptical array and in the intensity of the Friedel 
oscillations along the semimajor axis. For the 2DEG/achiral case, the calculations of the LDOS converged with fewer partial 
waves, i.e., the results were less sensitive to Imax under the same conditions as those used in the graphene calculations. This is 
most likely due to the fact that the /^^ partial wave is composed of both H^^"^ {kr) and H^^^ {kr), whereas the /^^ partial wave is 

composed of H^^\kr) in a 2DEG. It should also be mentioned that higher partial waves can lead to proximity resonances^^ as 
the scatterers are placed closer to one another; preliminary calculations indicate that these resonances, which are not observed 
for ^-wave scattering only [Imax = 0]. lead to large changes in the LDOS of graphene and should be experimentally observabl^. 
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Graphene case 



Achiral case 




FIG. 3: [Color Online] Theoretical calculations of fifai^.^l^ ^ ^ '^^^ ^ rectangular array of A/^ = 36 identical scatterers [V = 4 

eV, a = 10 A] for wave vectors k = 0.318 nm"^ [top, E = 0.21 eV for graphene], k = 0.419 nm"^ [middle, E = 0.2765 eV for graphene], 
and k = 0.485 nm~^ [bottom, E = 0.32 eV for graphene] for both a 2DEG [right, Eq. Jas}] and graphene [left, Eq. (|26|]. In all simulations, 

was set to zero within 35 A of each scatterer, si was calculated using Eq. ml in both graphene and the 2DEG, and Imax = 4 was 
P±K LJ 

used in all calculations in order to ensure a relative accuracy of calculations to within one percent. Dramatic differences in the quasiparticle 
interference patterns between the graphene [left] and 2DEG [right] cases are observed due to the interplay between multiple scattering and the 
resulting pseudospin rotations 



Finally, it should be noted that the scattering amplitudes in graphene, si in Eq. ([6]), are quite different than those found in 
2DEGs for scattering from a radial potential V(r), si^ac = ^^[j^^. These differences are attributable to the chiral nature of 

quasiparticles in graphene, where the wave function of both A and B lattice sites must be matched at each scatterer boundary. 
One consequence is that there exist conditions where l^-i | > |^o| even though ka <\, which is not true for si^ac- In these cases, 
taking into account only ^-wave scattering can lead to large errors in the calculated LDOS. This is illustrated for both an elliptical 
array [Fig. |4|A)] and a rectangular array [Fig. |4jB)], where the parameters used in the simulation are given in the caption of 
Fig. [4] In both cases, l^-i | ^ 7|^o|. which leads to large differences in the calculated LDOS if only 5'- wave scattering is considered; 
this is clearly illustrated in Figs. |4jA) and|4jB). 



V. CONCLUSIONS 



In this work, a theory for intravalley multiple scattering in single layer graphene was developed and used to calculate the 
total scattered wave function. Green's function, and the change in the local density of states (LDOS) in the presence of multiple 
scatterers. The scatterers were modeled by electrostatic step potentials that either raise or lower the potential within a radius a, 
with a being much larger than the C — C bond length, b = 1.42 A. Such scatterers could be experimentally realized by placing 
small metallic islands atop a graphene surface^ ^ In graphene, the shape, intensity, and pattern of quasiparticle interference 
were shown to be affected by the chiral nature of graphene compared with simple two-dimensional electron gases (2DEGs) or 2D 
achiral systems. The effects of these pseudospin rotations on the local density of states (LDOS) in graphene are similar to those 
predicted for the case of spin— orbit coupling in 2DEGS^^. Furthermore, the LDOS was found to be sensitive to the inclusion 
of higher partial waves if the distance between nearby scatterers was close to 2a, even in the ^--wave scattering scattering limit 
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FIG. 4: [Color Online] A comparison of the calculated fifai^^l^ [^^' ^'^^*] using either the first five partial waves [left, Imax = 4] or only the 

first or 5-wave scattering [Imax = 0. right] for [Fig. Qb)] a 200 A by 400 A rectangular array of A/^ = 36 identical scatterers [V = 4.93 eV, 
a= 10 A] and [Figs. |4[A)|4[C)] an elliptical array [semimajor and semiminor axis of 200 A and 100 A, respectively] composed of A/^ = 30 
identical scatterers [a = 10 Aand either (A) V = 4.93 eV or (C) V = 4 eV]. In all cases, fifai^^FN was set to zero within 35 A of the scatterers, 

and the choice of Imax = 4 [left] was found to provide a relative accuracy of one percent. For (A) and (B), the energies were chosen such 
that the p-wave scattering amplitude, l^i |, was larger than the 5-wave scattering amplitude, \so\, which gave (A) l^ol = 0.1429 < l^i | = 0.9975 
[k = 0.328 nm-\ E = 0.2163 eV] and (B) l^ol = 0.1237 < l^il = 0.9914 [k = 0.303 nm-^, E = 0.2 eV]. For the elliptical array in (C), 
however, even though |5o| = 0.9319 > = 0.0147 [k = 0.424 nm~^, E = 0.28 eV], there are still slight differences in near the 

array of scatterers between the calculations with Imax = 4 [left] and Imax = [right]. 



[i.e., I^ol ^ \^ly^o\ in Eq. The theoretical predictions in this work should be verifiable using STlVp^'^ and in atomicPl 

or microwave/optical'^^ experimental "simulations" of graphene systems. Finally, since the total Green's function is given by 
Eq. ([26]), the theory presented in this work could be used to calculate the effects of multiple scattering on the transport properties 
in graphene under conditions where intervalley scattering are negligible. 

The ultimate goal of this work was to lay out a formalism that enables the relatively fast calculation of the Green's function 
in single layer graphene for arbitrary scattering configurations and that elucidates the role that pseudospin interference has on 
the LDOS in graphene. A clear understanding of impurity scattering in graphene provides a step toward someday exploiting 
graphene 's unique properties to build new electronic devices. In the future, the theory in this work could be extended to scattering 
in multiple-layer graphene, nonzero bandgap graphene monolayers, higher energies, and to intervalley scattering. IMore realistic 
models for the scatterers could also be considered, where only the individual scattering amplitudes/phase shifts are needed as 
input in Eq. (26). Finally, resonances corresponding to those predicted for neutrino billiards might also be observabld^ for 
scatterers that locally induce a nonzero bandgap. 
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Appendix A: LDOS for 2DEG/achiral system 

The ^— matrix for scattering of higher partial waves by point scatterers in a 2DEG/2D achiral system has been previously 
derived by Herscll^. In this case, the total Green's function in the presence of N scatterers is given by: 

^ ^ ^ 2m^s^^^ ^ 



+ L £ (^+ [Go{7,rj,E)\ lL [Go{7j,7',E)\ [Go{7,rj,E)\ l'+ [G(ry,/,£)J J 



j=u=i 



-i^Hl,'\k\f-r'\) + l^sll^Hl,'\kpj)G{fjy,E) 



j=U=l ^ ^ ^ L J/ 



where Go{f,r,E) = -i^H^^\k\r -r'\), pj = \f-fjl e'^^J = ^(^ZzMMiy^ k = and m is the mass of the scattered 

particle. 

In this case, there are N{2lmax~^ 1) unknowns, which can be solved for self-consistently using the following Foldy— Lax 
equations [for / = to / = Imax and n = I to n = N]\ 



(A2) 



Defining the following ll^ax + 1x1 vectors for each scatterer y = 1 to j = N: 



( G{rj,f,E) \ 
L+[G(o-, ?',£)] 
L_[G(r,-, ?',£)] 
Ll[G{rj,f,E)] 
LL[G{fj,f,E)] 



TG {7j,7',E) = 



L^^--'[G(fj,r,E)] 
\L'r^[G{rj,r',E)] ) 



,TGq {rj,r) = -i—^ 
2n 



{-i)-^H^ll{kp])e-^''. 
{-iM\kp'j)e''''^ 
{-i)-^H^l{kp',)e-^'''. 



(-i)lma.H\^) {kp'.)e"'"''^^'j 



where pj = — rj \ and e^'^i = ''jH'^^'y) jj^g following N{2l,nax + 1 ) x 1 column vectors can then be constructed: 



' TG (ri,r') \ 
TG {72,7') 



\tG {7n,7')J 



,TacGo{7') = 



/ TGo (n,r') \ 
TGo {72,7') 



\TGo {7n,7')/ 



(A3) 



(A4) 



Eqs. ( A2 ) can be written compactly for n = 1 to « = and for / = to / = l,nca as: 



TacG{7') = 1-TTac T,,Go(/) 



(A5) 
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where 1 is the N{2ljnax + 1) x N{2lfnax + 1) identity matrix and TTac is a N{2ljnax + 1) x N{2lmax + 1) matrix given by: 

^ TTaciruh) TTaciruh) ... TTac{rurN) 



TTac = 



\ TTac{rNA) TTacirN.h) TTac{rN,h) 



where 



TTac{r„,rj) 



Sq L+ 



'Ho\krnj) 



s[^^H^'\kr,j) 



SJ)f2 



H^'hkrnj) 



Finally, G(f,r',E) in Eq. ( Al ) can be written compactly as: 







(A6) 



hnax 



Hl^\krnj) 



s\iH^'\kr„j) 



2> [ff^"(fer„,-)' 



G(r,r',£) = -/-^^('^(^Ir-r'l) +GGac(r)facG^^(r') 
In 

= -/^//W(fc|r-r'|) + ^ac(?)(l-TTac) 'TacGo±^(?'; 



''max ' 



H^o\kr,j) 
H^'\kr„j] 



njj 



s?ld'\kr,j) 



(A8) 



where GGac('^) is a 1 x N{2l,nax + 1) vector, GGac(?) 



GGac (f,ri) GGac (rji) ■■■ GGac (7,7^) 



where 



(A9) 



The 2DEG Green's function, G{r,r',E) in Eq. (A8), was used to calculate clean"/- l\ 

Fig.|2|andFig.|3| for the 2DEG/achiral 



cases. 



Appendix B: Equivalence between Foldy— Lax method and the point scatterer model 



In condensed matter systems, scatterers are often modeled as simple point scatterers. In this model, the scatterers are described 
by a sum of effecitve 5-potentials, V{f) = TJj=i o) where Vj, which has units of energy x area, is the potential associated 

with scatterer j. In this case, the total Green's function in the presence of N point scatterers, which we denote by G^^^ir^f' ^E), 
is given by the Lippmann-Sch winger equation: 

G^^^(r,/,£) = Go^^^(r,/,£) + y'dV'Go_^^(r,/',£)y(r")Gl^(r",/,£) 



G^^^(r,/,£) = Go_^^(r,/,£)+£y,Go_^^(r,r,-,£)G^^-(r,-,/,£) 



(Bl) 



As in the previous section, the various G^_^^{rj,r' ,E) can be solved for self-consistently by the following equations for 7 = 1 to 

Glg{rj,r',E) = G^^^girj,?' ,E) +VjG,^^,^{rj,rj,E)Gl.irj,r' ,E) + ^ y„Go_^^(ry,r„,£)G^^.(r„,r',£) (B2) 



Unlike the self-consistent Foldy— La^PSEU equations in Eq. (18), G^^ijj^r'^E) in Eq. (B2) is scaled by (1 — Gq ^^{f j j ^ E)) , 



where Gq ^^{f j j ^ E) is the "renormalized" Green's function of Gq ^^{f^f'^E) at r = = fj where the logarithmic singularity 
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has been removed, i.e., G^^^^ijj.rj.E) = + [ln(f ) +r]) 



1 
1 



where J is a renormaUzation length, and 



±K 



7= 0.5772 is the Euler-Mascheroni constant. The solution to Eqs. (B2) for the values of G^^{fj^r' ^E) for j = 1 to j = N and 
therefore the total Green's function, G^^^ir^r' ^E), can be written compactly using the Foldy— Lax eqatuion for 5'- wave scattering 

[Eq. (26) with 



- 0] but with the 5'- wave scattering amplitude for the scatterer, ^q^^ , replaced by 



"^0 - 



—I 



Ahvp 



1 



i^(in(f)+r) 



(B3) 



It is easy to verify that the point scatterer scattering amplitude, Sq^ in Eq. (B3), satisfies the unitarity condition, Re[i^Q^^] 



. Thus the point- scatterer model generates the same Green's function found using the Foldy— Lax equations for ^-wave 



scattering but with the ^-wave scattering amplitude given by Eq. (B3 ). 
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